// HybridMonteCarlo implementation
/////////////////////////////////////////////////////////////////
#include <iomanip>

# include "HybridMonteCarlo.h"
# include "Phasespace.h"

//////////////////////////////////////////////////////////////////
// Hybrid Monte Carlo : constructors/destructors
//////////////////////////////////////////////////////////////////

HybridMonteCarlo::HybridMonteCarlo(int NInput, int NStep,
                                   double LStep, Function *Func, double lProp)/*, double *mass)*/
                                   : StepRegulator(lProp)
{
  numDim = NInput;
  mTemperature = 1.0;

  StepNumber = NStep;
  StepLength = LStep;

//  NMin = 10.0;

  /*newA(m_mass,numDim);

  for(int i = 1; i<=numDim; i++)
  {
    m_mass[i] = mass[i];
  }*/

  InitialConfig.Initialise(NInput, Func);
  EvolvedConfig.Initialise(NInput, Func);

  StoredFunction = Func;

  StepRegulator.SetTarget(StepLength); 
}


//////////////////////////////////////////////////////////////////
// Hybrid Monte Carlo : functions
//////////////////////////////////////////////////////////////////

void HybridMonteCarlo::CreateSamples(double *Input, int NSample, int nMin)
{
  StepRegulator.ResetAll();

  InitialConfig.SetPosition(Input);
 

  for(int n=1; n<=NSample; n++)	
  {
	  HamiltonianPath();

    if(nMin == 0)
    {
      StoredFunction->Result.resetSecondSol();
    }
    else
    {
      StoredFunction->Result.resetSecondSol(n, nMin);
    }

	  Metropolis();
	}

  EvolvedConfig.CopyFinalPosition(Input);
}


void HybridMonteCarlo::CreateSamplesDiagnosis(double *Input, int NSample, const char filename[],
                                              const char filecoord[])
{
  double accept;
  nReject = 0;

  StepRegulator.ResetAll();

  ofstream outputcoords(filecoord);
  outputcoords.precision(6);
  outputcoords << setw(4) << 0 << setw(13);

  for(int l = 1; l<=numDim; l++)
  {
    outputcoords << Input[l] << setw(14);
  }
  outputcoords << endl;

  // setting up printing out diagnosis
  ofstream funcfile(filename);
  funcfile.precision(6);
  funcfile << "MC Step | Start Epot | End Epot |  Best Epot  |   delta_t   |   Func Eval   |   Distance"
    << endl << endl;


  InitialConfig.SetPosition(Input);

  double bestFuncValueSoFar; // used purely for outputting to funcfile

  for(int n=1; n<=NSample; n++)	
  {
    funcfile << setw(4) << n-1
      << setw(13) << InitialConfig.PotentialEnergy();
 
	  HamiltonianPath();

    // update Input to contain best solution so far for the 
    // purpose of calling StoredFunction->distanceFromGlobalMinimum(Input) below
    bestFuncValueSoFar = BestSoFar(Input);

  	funcfile << setw(13) << EvolvedConfig.PotentialEnergy() <<
      setw(13) << BestSoFar() << setw(13) << StepLength << 
      setw(14) << StoredFunction->nCall << setw(17) << StoredFunction->distanceFromGlobalMinimum(Input) << endl; 

	  accept = Metropolis();

    if(accept == 1)
    {
      outputcoords << setw(4) << n << setw(13);

      for(int k = 1; k<=numDim; k++)
      {
        outputcoords << InitialConfig.Position[k] << setw(14);
      }
      outputcoords << endl;
    }
	}

  EvolvedConfig.CopyFinalPosition(Input);

  funcfile.close();
  outputcoords.close();
}


void HybridMonteCarlo::HamiltonianPath(void)
{	
  // draw momentum vector from normal distribution with variance T
  InitialConfig.RandomiseMomentum(mTemperature);//, m_mass);

  EvolvedConfig = InitialConfig;


  EvolvedConfig.IncrementMomentum(StepLength/2.0);

  // StepNumber = Nmd
  for(int n=1;n<=StepNumber-1;n++)
	{
    EvolvedConfig.IncrementPosition(StepLength);//, m_mass);
    EvolvedConfig.PotentialEnergy();

	  EvolvedConfig.IncrementMomentum(StepLength);
	}

  EvolvedConfig.IncrementPosition(StepLength);//, m_mass);
  EvolvedConfig.PotentialEnergy();

  EvolvedConfig.IncrementMomentum(StepLength/2.0); 
}


int HybridMonteCarlo::Metropolis(void)
{
  double MetropFactor = InitialConfig.Energy() - EvolvedConfig.Energy();
  double diff = EvolvedConfig.PotentialEnergy() - StoredFunction->Result.getBestValue();

	MetropFactor /= mTemperature;

  if(exp(MetropFactor) >= ran1())	
  {
    //cout << *this;

    InitialConfig = EvolvedConfig;
    
    StepRegulator.GoodStep();

    return (1);
  }
  else
	{
	  //cout << "*" << endl;
	  StepRegulator.BadStep();
  	
    nReject++;

	  return (0);
	}
}


double HybridMonteCarlo::BestSoFar(double *solution)
{
  double x;

  StoredFunction->Result.Best(x, solution);

  return (x);
}


double HybridMonteCarlo::BestSoFar()
{
  double x;

  StoredFunction->Result.Best(x);

  return (x);
}


// special file to see what is going on for each trajectory
void HybridMonteCarlo::diagnoseTrajectory(string filename)
{

}

//////////////////////////////////////////////////////////////////
// HybridMonteCarlo : friends
//////////////////////////////////////////////////////////////////

ostream &operator<<(ostream &os, HybridMonteCarlo &X)
{
os << "(T=" << X.mTemperature << ")" << endl;
os << "Present state : ";

//os << "SE = " << X.InitialConfig.Energy() << "  ";
//os << "EE = " << X.EvolvedConfig.Energy() << "  ";
os << "PE = " << X.EvolvedConfig.PotentialEnergy() 
	<< "  Best = " << X.StoredFunction->Result.BestChiSq() << endl;

return ( os );
}
